1 Introduction

This webbook contains all the code used for data analysis in study of the individual-level metagenomic data of Podarcis muralis and Podarcis liolepis lizards from different environments during an experimental setup.

1.1 Prepare the R environment

1.1.1 Environment

To reproduce all the analyses locally, clone this repository in your computer using:

RStudio > New Project > Version Control > Git

And indicating the following git repository:

https://github.com/alberdilab/Podarcis_adaptation.git

Once the R project has been created, follow the instructions and code chunks shown in this webbook.

1.1.2 Libraries

The following R packages are required for the data analysis.

# Base
library(R.utils)
library(knitr)
library(tidyverse)
library(devtools)
library(tinytable)
library(janitor)

# For tree handling
library(ape)
library(phyloseq)
library(phytools)

# For plotting
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(ggnewscale)
library(gridExtra)
library(ggtreeExtra)
library(ggtree)
library(ggh4x)

# For statistics
library(spaa)
library(vegan)
library(Rtsne)
library(geiger)
library(hilldiv2)
library(distillR)
library(broom.mixed)
library(corrplot)
library(nlme)
library(pairwiseAdonis)

2 Prepare data

2.1 Load data

Load the original data files outputted by the bioinformatic pipeline.

2.1.1 Sample metadata

sample_metadata <- read_tsv("data/metadata.tsv")

2.1.2 Read counts

read_counts_raw <- read_tsv("data/genome.count.tsv") %>%
    rename(genome=1)

#Transformation of read_counts to combine data from both sequence rounds
merge_and_rename <- function(read_counts_raw) {
  read_counts_raw %>%
    # Gather the columns into long format
    pivot_longer(cols = -genome, names_to = "col") %>%
    # Extract prefix
    mutate(prefix = gsub("^(.*?)_.*", "\\1", col)) %>%
    # Group by prefix and genome, then summarize
    group_by(prefix, genome) %>%
    summarize(value = sum(value)) %>%
    # Spread the data back to wide format
    pivot_wider(names_from = prefix, values_from = value)
}

read_counts <- merge_and_rename(read_counts_raw)

2.1.3 Genome base hits

genome_hits_raw <- read_tsv("data/genome.covered_bases.tsv") %>%
    rename(genome=1)

#Transformation of genome_hits to combine data from both sequence rounds
merge_and_rename <- function(genome_hits_raw) {
  genome_hits_raw %>%
    # Gather the columns into long format
    pivot_longer(cols = -genome, names_to = "col") %>%
    # Extract prefix
    mutate(prefix = gsub("^(.*?)_.*", "\\1", col)) %>%
    # Group by prefix and genome, then summarize
    group_by(prefix, genome) %>%
    summarize(value = sum(value)) %>%
    # Spread the data back to wide format
    pivot_wider(names_from = prefix, values_from = value)
}

genome_hits <- merge_and_rename(genome_hits_raw)

2.1.4 Genome taxonomy

genome_taxonomy <- read_tsv("data/gtdbtk.summary.tsv") %>%
  select(mag_id = user_genome, classification) %>%
  separate(
    classification,
    into = c("domain", "phylum", "class", "order", "family", "genus", "species"),
    sep = ";") %>%
    rename(genome=1)

2.1.5 Genome quality

genome_quality <- read_tsv("data/quality_report.tsv") %>%
  select(
    genome = Name, 
    completeness = Completeness, 
    contamination = Contamination,
    length = Genome_Size, 
    gc = GC_Content
  )
genome_quality<-genome_quality %>% 
  mutate (genome=str_remove_all(genome,".fa"))

#Filter MAGs with over 70% completeness and less than 10% contamination
genome_quality <- genome_quality %>%
  filter(completeness > 70 & contamination < 10)

2.1.6 Genome tree

genome_tree <- read_tree("data/gtdbtk.backbone.bac120.classify.tree")
genome_tree$tip.label <- str_replace_all(genome_tree$tip.label,"'", "") #remove single quotes in MAG names

#Filter genome_taxonomy to keep MAGs with over 70% completeness and less than 10% contamination
genome_taxonomy <- genome_taxonomy %>%
  semi_join(genome_quality, by = "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_taxonomy$genome) # keep only MAG tips

2.1.7 Genome annotations

genome_annotations <- read_tsv("data/annotations.tsv.xz") %>%
    rename(gene=1, genome=2, contig=3)

#Filter only the MAGs with over 70% completeness and less than 10% contamination
genome_annotations <- genome_annotations %>%
  semi_join(genome_quality, by = "genome")

2.2 Create working objects

Transform the original data files into working objects for downstream analyses.

2.2.1 Merge genome taxonomy and quality

genome_metadata <- genome_taxonomy %>%
  inner_join(genome_quality,by=join_by(genome==genome)) #join quality

2.2.2 Calculate genome coverage

#Filter genome_hits for the MAGs with over 70% completeness and less than 10% contamination
genome_hits <- genome_hits %>%
  semi_join(genome_quality, by = "genome")

genome_coverage <- genome_hits %>%
  mutate(across(where(is.numeric), ~ ./genome_metadata$length))

2.3 Filtering

4 samples are removed from the analysis due to their low sequencing depth.

#Counts_raw
columns_to_exclude <- c("AD16","AD23", "AD25", "AD91") # Columns to exclude
read_counts <- read_counts %>%
  select(-columns_to_exclude)

#Metadata
sample_metadata <- sample_metadata %>%
  filter(Tube_code != "AD16") %>%
  filter(Tube_code != "AD23") %>%
  filter(Tube_code != "AD25") %>%
  filter(Tube_code != "AD91")

#Coverage_table
columns_to_exclude <- c("AD16", "AD23", "AD25", "AD91")  # Columns to exclude
genome_coverage <- genome_coverage %>%
  select(-columns_to_exclude)

##Removal of samples of individual LI1_2nd_6
#Counts_raw
columns_to_exclude <- c("AFV13","AD06", "AD31", "AE03", "AF12") # Columns to exclude
read_counts <- read_counts %>%
  select(-columns_to_exclude)

#Metadata
sample_metadata <- sample_metadata %>%
  filter(Tube_code != "AFV13") %>%
  filter(Tube_code != "AD06") %>%
  filter(Tube_code != "AD31") %>%
  filter(Tube_code != "AE03") %>%
  filter(Tube_code != "AF12") 

#Coverage_table
columns_to_exclude <- c("AFV13","AD06", "AD31", "AE03", "AF12")  # Columns to exclude
genome_coverage <- genome_coverage %>%
  select(-columns_to_exclude)

2.3.1 Filter reads by coverage

#Filter read_counts for the MAGs with over 70% completeness and less than 10% contamination
read_counts <- read_counts %>%
  semi_join(genome_quality, by = "genome")

min_coverage=0.3
read_counts_filt <- genome_coverage %>%
  mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
  mutate(across(-1, ~ . * read_counts[[cur_column()]])) 

2.3.2 Transform reads into genome counts

readlength=150
genome_counts <- read_counts %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

readlength=150
genome_counts_filt <- read_counts_filt %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

2.3.3 Distill annotations into GIFTs

genome_gifts <- distill(genome_annotations,GIFT_db,genomecol=2,annotcol=c(9,10,19))

2.4 Load data statistics

2.4.1 Raw reads

raw_reads <-
  "data/multiqc_general_stats_2.txt" %>%
  read_tsv() %>%
  select(
    sample = Sample,
    raw_reads = `total_sequences`
  ) %>%
  mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
  group_by(sample_prefix) %>%
  summarize(raw_reads = sum(raw_reads, na.rm = TRUE)) %>%
  rename(sample = sample_prefix)  %>%
  mutate(sample = str_remove(sample, "^fastp \\|\\s*"))

2.4.2 Quality-filtered reads

fastp_reads <-
  "data/multiqc_general_stats.txt" %>%
  read_tsv() %>%
  select(
    sample = Sample,
    trimmed_reads = `total_sequences`
  ) %>%
  mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
  group_by(sample_prefix) %>%
  summarize(trimmed_reads = sum(trimmed_reads, na.rm = TRUE)) %>%
  rename(sample = sample_prefix) %>% 
  filter(!str_detect(sample, "nonlizard \\|")) %>%
  filter(!str_detect(sample, "lizard \\|")) %>%
  filter(!str_detect(sample, "refseq500 \\|"))  %>%
  mutate(sample = str_remove(sample, "^fastp \\|\\s*"))

2.4.3 Host-mapped reads

host_mapped <-
  "data/multiqc_general_stats.txt" %>%
  read_tsv() %>%
  filter(str_detect(Sample, "lizard")) %>%
  select(
    sample = Sample,
    host_mapped = `reads_mapped`,
    mapping_total = `raw_total_sequences`
  ) %>%
  mutate(
    host_unmapped = mapping_total - host_mapped
  ) %>%
  filter(!is.na(host_mapped)) %>%
  separate(
    col = sample,
    into = c("host_name", "sample"),
    sep = " \\| "
  ) %>%
  rename(mapped = host_mapped, unmapped = host_unmapped) %>%
  select(-mapping_total) %>%
  pivot_longer(-host_name:-sample) %>%
  mutate(
    name = str_glue("{name}_{host_name}")
  ) %>%
  select(-host_name) %>%
  pivot_wider() %>%
  mutate(sample = str_extract(sample, "^[^_]+")) %>%
  group_by(sample) %>%
  summarize(
    mapped_lizard = sum(mapped_lizard),
    unmapped_lizard = sum(unmapped_lizard)
  ) 

2.4.4 Prokaryotic fraction

singlem <-
  "data/singleM.tsv" %>%
  read_tsv() %>%
  distinct() %>%
  mutate(
    sample = sample %>% str_remove_all("_1$"),
   read_fraction = read_fraction %>% str_remove("%") %>% as.numeric(),
   read_fraction = read_fraction / 100
  ) %>%
  select(
   sample,
   singlem_prokaryotic_bases = bacterial_archaeal_bases,
   singlem_metagenome_size = metagenome_size,
   singlem_read_fraction = read_fraction,
  ) %>%
mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
  group_by(sample_prefix) %>%
  summarize(
    singlem_prokaryotic_bases = sum(singlem_prokaryotic_bases),
    singlem_metagenome_size = sum(singlem_metagenome_size),
    singlem_read_fraction = mean(singlem_read_fraction)
  ) %>%
  rename(sample = sample_prefix)

2.4.5 MAG-mapped reads

mag_mapping <-
  read_tsv("data/contig.count.tsv.xz") %>%
  pivot_longer(-sequence_id) %>%
  summarise(value = sum(value), .by = "name") %>%
  rename(sample = name, mapped_mags = value) %>%
  mutate(sample_prefix = str_extract(sample, "^[^_]+")) %>%
    group_by(sample_prefix) %>%
    summarize(
      mapped_mags = sum(mapped_mags)
    ) %>%
    rename(sample = sample_prefix)

2.4.6 Wrap data statistics

data_stats <- raw_reads %>%
  left_join(fastp_reads) %>%
  left_join(host_mapped) %>%
  left_join(singlem) %>%
  left_join(mag_mapping)

data_stats<- data_stats %>%
  filter(!str_detect(sample, "nonlizard \\|")) %>%
  filter(!str_detect(sample, "lizard \\|")) %>%
  filter(!str_detect(sample, "refseq500 \\|"))

2.5 Prepare color scheme

AlberdiLab projects use unified color schemes developed for the Earth Hologenome Initiative, to facilitate figure interpretation.

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    arrange(match(genome, genome_tree$tip.label)) %>%
    select(phylum, colors) %>% 
    unique() %>%
    arrange(phylum) %>%
    pull(colors, name=phylum)

2.6 Wrap working objects

All working objects are wrapped into a single Rdata object to facilitate downstream usage.

save(sample_metadata, 
     genome_metadata, 
     read_counts, 
     genome_counts, 
     genome_counts_filt, 
     genome_tree,
     genome_gifts, 
     phylum_colors,
     data_stats,
     file = "data/data.Rdata")

3 Data statistics

load("data/data.Rdata")

3.1 Sequencing reads statistics

data_stats$raw_reads %>% sum()

[1] 4995505290

data_stats$raw_reads %>% mean()

[1] 28875753

data_stats$raw_reads %>% sd()

[1] 13447526

3.2 DNA fractions

#Overall
data_stats %>%
    mutate(mapped_perc=mapped_mags/trimmed_reads) %>%
    summarise(mean=mean(mapped_perc),sd=sd(mapped_perc))

# A tibble: 1 × 2
   mean    sd
  <dbl> <dbl>
1 0.438 0.207

data_stats %>%
  mutate(
    low_quality = raw_reads - trimmed_reads,
    unmapped_reads = trimmed_reads - mapped_lizard - mapped_mags
  ) %>%
  select(sample, low_quality, mapped_lizard, mapped_mags, unmapped_reads) %>%
  pivot_longer(-sample) %>%
  mutate(name=factor(name,levels=c("low_quality","mapped_lizard","unmapped_reads","mapped_mags"))) %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(!is.na(time_point)) %>%
  ggplot(aes(x = sample, y = value, fill = name)) +
      geom_bar(stat = "identity", position = "fill") +
      scale_fill_manual(name="Sequence type",
                    breaks=c("low_quality","mapped_lizard","unmapped_reads","mapped_mags"),
                    labels=c("Low quality","Mapped to host","Unmapped","Mapped to MAGs"),
                    values=c("#CCCCCC", "#bcdee1", "#d8b8a3","#93655c"))+
      facet_grid(~time_point, scales = "free", labeller = labeller(time_point=c("0_Wild"="Wild", "1_Acclimation"="Acclimation", "2_Antibiotics"="Antibiotics",
                                                                                "3_Transplant1"="Transplant1", "4_Transplant2"="Transplant2", "5_Post-FMT1"="Post1",
                                                                                "6_Post-FMT2"="Post2"))) +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 90, hjust = 1, size=0)) +
      labs(y="DNA sequence fraction",x="Samples")

3.3 Recovered microbial fraction

data_stats %>%
  mutate(
    unmapped_reads = trimmed_reads - mapped_lizard - mapped_mags,
    mag_proportion = mapped_mags / (mapped_mags + unmapped_reads),
    singlem_read_fraction = singlem_read_fraction
  ) %>%
  select(sample, mag_proportion, singlem_read_fraction) %>%
  mutate(
    mag_proportion = if_else(singlem_read_fraction == 0, 0, mag_proportion),
    singlem_read_fraction = if_else(singlem_read_fraction == 0, NA, singlem_read_fraction),
    singlem_read_fraction = if_else(singlem_read_fraction < mag_proportion, NA, singlem_read_fraction),
    singlem_read_fraction = if_else(singlem_read_fraction > 1, 1, singlem_read_fraction)
  ) %>%
  pivot_longer(-sample, names_to = "proportion", values_to = "value") %>%
  mutate(
    proportion = factor(
      proportion,
      levels = c("mag_proportion", "singlem_read_fraction")
    )
  ) %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(!is.na(time_point)) %>%
  ggplot(aes(x = sample, y = value, color = proportion)) +
      geom_line(aes(group = sample), color = "#f8a538") +
      geom_point() +
      scale_color_manual(name="Proportion",
                    breaks=c("mag_proportion","singlem_read_fraction"),
                    labels=c("Recovered","Estimated"),
                    values=c("#52e1e8", "#876b53"))+
      theme_minimal() +
      facet_grid(~time_point, scales = "free", space="free", labeller = labeller(time_point=c("0_Wild"="Wild", "1_Acclimation"="Acclimation", "2_Antibiotics"="Antibiotics",
                                                                                "3_Transplant1"="Transplant1", "4_Transplant2"="Transplant2", "5_Post-FMT1"="Post1",
                                                                                "6_Post-FMT2"="Post2")))+
      labs(y = "Samples", x = "Prokaryotic fraction") +
      scale_y_continuous(limits = c(0, 1)) +
      theme(
        axis.text.y = element_text(size = 4),
        axis.text.x = element_text( angle = 90, vjust = 0.5, hjust = 1, size = 0),
        legend.position = "right"
      )

3.3.1 Domain-adjusted mapping rate (DAMR)

data_stats %>%
  mutate(
    unmapped_reads = trimmed_reads - mapped_lizard - mapped_mags,
    mag_proportion = mapped_mags / (mapped_mags + unmapped_reads),
    singlem_read_fraction = singlem_read_fraction
  ) %>%
  mutate(damr=pmin(1, mag_proportion/singlem_read_fraction)) %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(!is.na(time_point)) %>%
  select(sample,damr, time_point, species, type) %>%
  #group_by(type) %>%
  #summarise(mean=mean(damr),sd=sd(damr)) %>%
  tt()

sample	damr	time_point	species	type
AC79	1.0000000	1_Acclimation	Podarcis_muralis	Control
AC80	1.0000000	1_Acclimation	Podarcis_muralis	Treatment
AC81	0.9813013	1_Acclimation	Podarcis_muralis	Treatment
AC82	0.7490771	1_Acclimation	Podarcis_muralis	Treatment
AC83	1.0000000	1_Acclimation	Podarcis_muralis	Control
AC84	1.0000000	1_Acclimation	Podarcis_muralis	Control
AC85	0.8363777	1_Acclimation	Podarcis_muralis	Control
AC86	0.8800275	1_Acclimation	Podarcis_muralis	Control
AC87	0.9672854	1_Acclimation	Podarcis_muralis	Treatment
AC88	0.3142108	1_Acclimation	Podarcis_muralis	Control
AC89	0.9338962	1_Acclimation	Podarcis_muralis	Treatment
AC90	1.0000000	1_Acclimation	Podarcis_muralis	Control
AC91	1.0000000	1_Acclimation	Podarcis_muralis	Control
AC92	1.0000000	1_Acclimation	Podarcis_muralis	Control
AC93	1.0000000	1_Acclimation	Podarcis_muralis	Treatment
AC94	0.9159348	1_Acclimation	Podarcis_muralis	Treatment
AC95	0.9239904	1_Acclimation	Podarcis_muralis	Treatment
AC96	0.9425962	1_Acclimation	Podarcis_muralis	Treatment
AC97	0.9394886	1_Acclimation	Podarcis_liolepis	Hot_control
AC98	1.0000000	1_Acclimation	Podarcis_liolepis	Hot_control
AC99	1.0000000	1_Acclimation	Podarcis_liolepis	Hot_control
AD01	1.0000000	1_Acclimation	Podarcis_liolepis	Hot_control
AD02	1.0000000	1_Acclimation	Podarcis_liolepis	Hot_control
AD03	1.0000000	1_Acclimation	Podarcis_liolepis	Hot_control
AD04	0.8896057	1_Acclimation	Podarcis_liolepis	Hot_control
AD05	1.0000000	1_Acclimation	Podarcis_liolepis	Hot_control
AD07	0.9515379	1_Acclimation	Podarcis_liolepis	Hot_control
AD08	0.7045003	2_Antibiotics	Podarcis_muralis	Control
AD09	1.0000000	2_Antibiotics	Podarcis_muralis	Control
AD10	0.9598720	2_Antibiotics	Podarcis_muralis	Control
AD11	0.6230176	2_Antibiotics	Podarcis_muralis	Control
AD12	0.7487463	2_Antibiotics	Podarcis_muralis	Control
AD13	1.0000000	2_Antibiotics	Podarcis_muralis	Control
AD14	0.4536529	2_Antibiotics	Podarcis_muralis	Control
AD15	1.0000000	2_Antibiotics	Podarcis_muralis	Control
AD17	0.8774871	2_Antibiotics	Podarcis_muralis	Treatment
AD18	0.4991585	2_Antibiotics	Podarcis_muralis	Treatment
AD19	0.4060523	2_Antibiotics	Podarcis_muralis	Treatment
AD20	0.5227482	2_Antibiotics	Podarcis_muralis	Treatment
AD21	0.6551882	2_Antibiotics	Podarcis_muralis	Treatment
AD22	0.9106725	2_Antibiotics	Podarcis_muralis	Treatment
AD24	0.5972915	2_Antibiotics	Podarcis_muralis	Treatment
AD26	0.9204825	2_Antibiotics	Podarcis_liolepis	Hot_control
AD27	0.5215059	2_Antibiotics	Podarcis_liolepis	Hot_control
AD28	0.8180607	2_Antibiotics	Podarcis_liolepis	Hot_control
AD29	0.7175978	2_Antibiotics	Podarcis_liolepis	Hot_control
AD30	0.9937574	2_Antibiotics	Podarcis_liolepis	Hot_control
AD32	0.6710817	2_Antibiotics	Podarcis_liolepis	Hot_control
AD33	0.6778368	2_Antibiotics	Podarcis_liolepis	Hot_control
AD34	0.4642326	2_Antibiotics	Podarcis_liolepis	Hot_control
AD36	0.8288365	3_Transplant1	Podarcis_muralis	Control
AD37	0.9530848	3_Transplant1	Podarcis_muralis	Control
AD38	0.9247925	3_Transplant1	Podarcis_muralis	Control
AD39	1.0000000	3_Transplant1	Podarcis_muralis	Control
AD40	0.9991552	3_Transplant1	Podarcis_muralis	Control
AD41	0.9550964	3_Transplant1	Podarcis_muralis	Control
AD42	1.0000000	3_Transplant1	Podarcis_muralis	Control
AD43	1.0000000	3_Transplant1	Podarcis_muralis	Control
AD44	0.9108807	3_Transplant1	Podarcis_muralis	Control
AD46	0.8708468	3_Transplant1	Podarcis_muralis	Treatment
AD47	0.8328653	3_Transplant1	Podarcis_muralis	Treatment
AD49	0.8940294	3_Transplant1	Podarcis_muralis	Treatment
AD50	0.8757711	3_Transplant1	Podarcis_muralis	Treatment
AD51	0.8479193	3_Transplant1	Podarcis_muralis	Treatment
AD52	0.8596870	3_Transplant1	Podarcis_muralis	Treatment
AD53	0.8528974	3_Transplant1	Podarcis_muralis	Treatment
AD54	1.0000000	3_Transplant1	Podarcis_liolepis	Hot_control
AD55	1.0000000	4_Transplant2	Podarcis_muralis	Control
AD56	0.9327538	4_Transplant2	Podarcis_muralis	Control
AD57	1.0000000	4_Transplant2	Podarcis_muralis	Control
AD58	1.0000000	4_Transplant2	Podarcis_muralis	Control
AD59	1.0000000	4_Transplant2	Podarcis_muralis	Control
AD60	0.7988844	4_Transplant2	Podarcis_muralis	Control
AD61	1.0000000	4_Transplant2	Podarcis_muralis	Control
AD62	1.0000000	4_Transplant2	Podarcis_muralis	Control
AD63	0.8992048	4_Transplant2	Podarcis_muralis	Control
AD65	1.0000000	4_Transplant2	Podarcis_muralis	Treatment
AD66	0.8821027	4_Transplant2	Podarcis_muralis	Treatment
AD68	0.8820274	4_Transplant2	Podarcis_muralis	Treatment
AD70	0.9590793	4_Transplant2	Podarcis_muralis	Treatment
AD71	0.8360540	4_Transplant2	Podarcis_muralis	Treatment
AD72	0.8835253	4_Transplant2	Podarcis_muralis	Treatment
AD73	1.0000000	4_Transplant2	Podarcis_liolepis	Hot_control
AD74	0.9017271	5_Post-FMT1	Podarcis_muralis	Control
AD75	1.0000000	5_Post-FMT1	Podarcis_muralis	Treatment
AD76	0.9862581	5_Post-FMT1	Podarcis_muralis	Treatment
AD77	0.8963376	5_Post-FMT1	Podarcis_muralis	Treatment
AD78	1.0000000	5_Post-FMT1	Podarcis_muralis	Treatment
AD79	1.0000000	5_Post-FMT1	Podarcis_muralis	Control
AD80	0.8808155	5_Post-FMT1	Podarcis_muralis	Treatment
AD81	1.0000000	5_Post-FMT1	Podarcis_muralis	Control
AD82	1.0000000	5_Post-FMT1	Podarcis_muralis	Control
AD83	1.0000000	5_Post-FMT1	Podarcis_muralis	Treatment
AD84	1.0000000	5_Post-FMT1	Podarcis_muralis	Treatment
AD85	1.0000000	5_Post-FMT1	Podarcis_muralis	Control
AD86	0.9889084	5_Post-FMT1	Podarcis_muralis	Control
AD87	1.0000000	5_Post-FMT1	Podarcis_muralis	Control
AD88	1.0000000	5_Post-FMT1	Podarcis_muralis	Treatment
AD89	0.9562489	5_Post-FMT1	Podarcis_muralis	Control
AD90	1.0000000	5_Post-FMT1	Podarcis_muralis	Control
AD93	1.0000000	5_Post-FMT1	Podarcis_liolepis	Hot_control
AD94	0.7947750	5_Post-FMT1	Podarcis_liolepis	Hot_control
AD95	0.9506977	5_Post-FMT1	Podarcis_liolepis	Hot_control
AD96	0.9171413	5_Post-FMT1	Podarcis_liolepis	Hot_control
AD97	0.9257685	5_Post-FMT1	Podarcis_liolepis	Hot_control
AD98	0.9048349	5_Post-FMT1	Podarcis_liolepis	Hot_control
AD99	0.9966152	5_Post-FMT1	Podarcis_liolepis	Hot_control
AE01	1.0000000	5_Post-FMT1	Podarcis_liolepis	Hot_control
AE02	1.0000000	5_Post-FMT1	Podarcis_liolepis	Hot_control
AE04	1.0000000	6_Post-FMT2	Podarcis_muralis	Treatment
AE05	1.0000000	6_Post-FMT2	Podarcis_muralis	Treatment
AE06	1.0000000	6_Post-FMT2	Podarcis_muralis	Control
AE07	0.7797794	6_Post-FMT2	Podarcis_muralis	Control
AE08	1.0000000	6_Post-FMT2	Podarcis_muralis	Control
AE09	1.0000000	6_Post-FMT2	Podarcis_muralis	Control
AE91	0.9259146	6_Post-FMT2	Podarcis_muralis	Treatment
AE92	1.0000000	6_Post-FMT2	Podarcis_muralis	Control
AE93	1.0000000	6_Post-FMT2	Podarcis_muralis	Control
AE94	0.9152960	6_Post-FMT2	Podarcis_muralis	Treatment
AE95	1.0000000	6_Post-FMT2	Podarcis_muralis	Treatment
AE96	0.7786112	6_Post-FMT2	Podarcis_muralis	Treatment
AE97	0.9365037	6_Post-FMT2	Podarcis_muralis	Control
AE98	1.0000000	6_Post-FMT2	Podarcis_muralis	Control
AE99	0.8513354	6_Post-FMT2	Podarcis_muralis	Treatment
AF01	0.9046369	6_Post-FMT2	Podarcis_muralis	Treatment
AF02	1.0000000	6_Post-FMT2	Podarcis_muralis	Control
AF03	0.9030087	6_Post-FMT2	Podarcis_muralis	Treatment
AF04	1.0000000	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF05	1.0000000	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF06	1.0000000	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF07	0.9526443	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF08	1.0000000	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF09	1.0000000	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF10	0.9716162	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF11	0.9654152	6_Post-FMT2	Podarcis_liolepis	Hot_control
AF13	0.8818742	6_Post-FMT2	Podarcis_liolepis	Hot_control
AFU87	0.8810513	0_Wild	Podarcis_muralis	Control
AFU88	0.8317800	0_Wild	Podarcis_muralis	Treatment
AFU91	0.8923699	0_Wild	Podarcis_muralis	Treatment
AFU92	0.8094776	0_Wild	Podarcis_muralis	Treatment
AFU93	0.8517368	0_Wild	Podarcis_muralis	Control
AFU94	0.8325385	0_Wild	Podarcis_muralis	Treatment
AFU95	0.8419270	0_Wild	Podarcis_muralis	Treatment
AFU96	0.8326820	0_Wild	Podarcis_muralis	Control
AFU97	0.8107271	0_Wild	Podarcis_muralis	Treatment
AFU98	0.7506522	0_Wild	Podarcis_muralis	Control
AFU99	0.8582371	0_Wild	Podarcis_muralis	Treatment
AFV01	0.9331539	0_Wild	Podarcis_muralis	Treatment
AFV02	0.8316460	0_Wild	Podarcis_muralis	Treatment
AFV03	0.8752591	0_Wild	Podarcis_muralis	Control
AFV04	0.9180285	0_Wild	Podarcis_muralis	Control
AFV05	1.0000000	0_Wild	Podarcis_muralis	Control
AFV06	1.0000000	0_Wild	Podarcis_muralis	Control
AFV07	0.8460805	0_Wild	Podarcis_muralis	Control
AFV08	0.7497043	0_Wild	Podarcis_liolepis	Hot_control
AFV09	0.5412999	0_Wild	Podarcis_liolepis	Hot_control
AFV10	0.8002499	0_Wild	Podarcis_liolepis	Hot_control
AFV11	0.8225298	0_Wild	Podarcis_liolepis	Hot_control
AFV12	0.7925988	0_Wild	Podarcis_liolepis	Hot_control
AFV14	0.8106269	0_Wild	Podarcis_liolepis	Hot_control
AFV15	0.9691106	0_Wild	Podarcis_liolepis	Hot_control
AFV16	0.8218990	0_Wild	Podarcis_liolepis	Hot_control
AFV17	0.8091152	0_Wild	Podarcis_liolepis	Hot_control

4 MAG catalogue

load("data/data.Rdata")

4.1 Genome phylogeny

# Generate the phylum color heatmap
phylum_heatmap <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    arrange(match(genome, genome_tree$tip.label)) %>%
    select(genome,phylum) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome")

# Generate  basal tree
circular_tree <- force.ultrametric(genome_tree, method="extend") %>% # extend to ultrametric for the sake of visualisation
    ggtree(., layout="fan", open.angle=10, size=0.5)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

# Add phylum ring
circular_tree <- gheatmap(circular_tree, phylum_heatmap, offset=0.55, width=0.1, colnames=FALSE) +
        scale_fill_manual(values=phylum_colors) +
        geom_tiplab2(size=1, hjust=-0.1) +
        theme(legend.position = "none", plot.margin = margin(0, 0, 0, 0), panel.margin = margin(0, 0, 0, 0))

# Flush color scale to enable a new color scheme in the next ring
circular_tree <- circular_tree + new_scale_fill()

# Add completeness ring
circular_tree <- circular_tree +
        new_scale_fill() +
        scale_fill_gradient(low = "#d1f4ba", high = "#f4baba") +
        geom_fruit(
                data=genome_metadata,
                geom=geom_bar,
                mapping = aes(x=completeness, y=genome, fill=contamination),
                offset = 0.55,
                orientation="y",
              stat="identity")

# Add genome-size ring
circular_tree <-  circular_tree +
        new_scale_fill() +
        scale_fill_manual(values = "#cccccc") +
        geom_fruit(
             data=genome_metadata,
             geom=geom_bar,
             mapping = aes(x=length, y=genome),
                 offset = 0.05,
                 orientation="y",
         stat="identity")

# Add text
circular_tree <-  circular_tree +
        annotate('text', x=3.4, y=0, label='              Phylum', family='ArialMT', size=3.5) +
        annotate('text', x=4.2, y=0, label='                         Genome quality', family='ArialMT', size=3.5) +
        annotate('text', x=4.7, y=0, label='                      Genome size', family='ArialMT', size=3.5)

#Plot circular tree
circular_tree %>% open_tree(30) %>% rotate_tree(90)

4.2 Genome quality

genome_metadata$completeness %>% mean()

[1] 92.98658

genome_metadata$completeness %>% sd()

[1] 7.170467

genome_metadata$contamination %>% mean()

[1] 2.154888

genome_metadata$contamination %>% sd()

[1] 2.194372

#Generate quality biplot
genome_biplot <- genome_metadata %>%
  select(c(genome,domain,phylum,completeness,contamination,length)) %>%
  arrange(match(genome, rev(genome_tree$tip.label))) %>% #sort MAGs according to phylogenetic tree
  ggplot(aes(x=completeness,y=contamination,size=length,color=phylum)) +
              geom_point(alpha=0.7) +
                    ylim(c(10,0)) +
                    scale_color_manual(values=phylum_colors) +
                    labs(y= "Contamination", x = "Completeness") +
                    theme_classic() +
                    theme(legend.position = "none")

#Generate contamination boxplot
genome_contamination <- genome_metadata %>%
            ggplot(aes(y=contamination)) +
                    ylim(c(10,0)) +
                    geom_boxplot(colour = "#999999", fill="#cccccc") +
                    theme_void() +
                    theme(legend.position = "none",
                        axis.title.x = element_blank(),
                        axis.title.y = element_blank(),
                        axis.text.y=element_blank(),
                        axis.ticks.y=element_blank(),
                        axis.text.x=element_blank(),
                        axis.ticks.x=element_blank(),
                        plot.margin = unit(c(0, 0, 0.40, 0),"inches")) #add bottom-margin (top, right, bottom, left)

#Generate completeness boxplot
genome_completeness <- genome_metadata %>%
        ggplot(aes(x=completeness)) +
                xlim(c(50,100)) +
                geom_boxplot(colour = "#999999", fill="#cccccc") +
                theme_void() +
                theme(legend.position = "none",
                    axis.title.x = element_blank(),
                    axis.title.y = element_blank(),
                    axis.text.y=element_blank(),
                    axis.ticks.y=element_blank(),
                    axis.text.x=element_blank(),
                    axis.ticks.x=element_blank(),
                    plot.margin = unit(c(0, 0, 0, 0.50),"inches")) #add left-margin (top, right, bottom, left)

#Render composite figure
grid.arrange(grobs = list(genome_completeness,genome_biplot,genome_contamination),
        layout_matrix = rbind(c(1,1,1,1,1,1,1,1,1,1,1,4),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3),
                              c(2,2,2,2,2,2,2,2,2,2,2,3)))

4.3 Functional overview

# Aggregate basal GIFT into elements
function_table <- genome_gifts %>%
    to.elements(., GIFT_db)

# Generate  basal tree
function_tree <- force.ultrametric(genome_tree, method="extend") %>%
                ggtree(., size = 0.3) 

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#Add phylum colors next to the tree tips
function_tree <- gheatmap(function_tree, phylum_heatmap, offset=0, width=0.1, colnames=FALSE) +
            scale_fill_manual(values=phylum_colors) +
            labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
function_tree <- function_tree + new_scale_fill()

#Add functions heatmap
function_tree <- gheatmap(function_tree, function_table, offset=0.5, width=3.5, colnames=FALSE) +
            vexpand(.08) +
            coord_cartesian(clip = "off") +
            scale_fill_gradient(low = "#f4f4f4", high = "steelblue", na.value="white") +
            labs(fill="GIFT")

#Reset fill scale to use a different colour profile in the heatmap
function_tree <- function_tree + new_scale_fill()

# Add completeness barplots
function_tree <- function_tree +
            geom_fruit(data=genome_metadata,
            geom=geom_bar,
            grid.params=list(axis="x", text.size=2, nbreak = 1),
            axis.params=list(vline=TRUE),
            mapping = aes(x=length, y=genome, fill=completeness),
                 offset = 3.8,
                 orientation="y",
                 stat="identity") +
            scale_fill_gradient(low = "#cf8888", high = "#a2cc87") +
            labs(fill="Genome\ncompleteness")

function_tree

4.4 Functional ordination

# Generate the tSNE ordination
tSNE_function <- Rtsne(X=function_table, dims = 2, check_duplicates = FALSE)

# Plot the ordination
function_ordination <- tSNE_function$Y %>%
                as.data.frame() %>%
                mutate(genome=rownames(function_table)) %>%
                inner_join(genome_metadata, by="genome") %>%
                rename(tSNE1="V1", tSNE2="V2") %>%
                select(genome,phylum,tSNE1,tSNE2, length) %>%
                ggplot(aes(x = tSNE1, y = tSNE2, color = phylum, size=length))+
                            geom_point(shape=16, alpha=0.7) +
                            scale_color_manual(values=phylum_colors) +
                            theme_minimal() +
                labs(color="Phylum", size="Genome size") +
                guides(color = guide_legend(override.aes = list(size = 5))) # enlarge Phylum dots in legend

function_ordination

5 Community composition

load("data/data.Rdata")

5.1 Taxonomy overview

5.1.1 Stacked barplot

# Merge data frames based on sample
transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

merged_data<-genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., transplants_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  filter(count > 0) #filter 0 counts

ggplot(merged_data, aes(x=newID,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ time_point + type ,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    labs(fill="Phylum",y = "Relative abundance",x="Sample")+
    theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size=0))

5.1.1.1 Wild samples

merged_data  %>%
  filter(time_point=="0_Wild")  %>%
  ggplot(aes(x=newID,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ Population,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    labs(fill="Phylum",y = "Relative abundance",x="Sample")+
    theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size=0),
    strip.text.x = element_text(size = 12)
    )

5.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(total_mean=mean(relabun*100, na.rm=T),
              total_sd=sd(relabun*100, na.rm=T))  %>%
    mutate(total=str_c(round(total_mean,2),"±",round(total_sd,2))) %>% 
    arrange(-total_mean) %>% 
    dplyr::select(phylum,total) %>% 
    tt()

phylum	total
p__Bacteroidota	38.75±19.94
p__Bacillota_A	24.77±15.73
p__Bacillota	11.96±14.79
p__Pseudomonadota	9.41±15.86
p__Campylobacterota	5.49±9.42
p__Verrucomicrobiota	2.78±6.7
p__Desulfobacterota	2.37±3.68
p__Chlamydiota	1.1±6.08
p__Fusobacteriota	1.06±2.86
p__Cyanobacteriota	0.93±1.66
p__Bacillota_C	0.48±0.67
p__Spirochaetota	0.41±1.25
p__Bacillota_B	0.26±0.5
p__Actinomycetota	0.13±0.65
p__Elusimicrobiota	0.11±0.63

5.1.2.1 Wild samples

#Cold and wet population
phylum_summary %>%
    left_join(sample_metadata, by = join_by(sample == Tube_code))  %>%
    group_by(phylum, Population) %>%
    filter(Population=="Cold_wet") %>%
    filter(time_point=="0_Wild") %>%
    summarise(total_mean=mean(relabun*100, na.rm=T),
              total_sd=sd(relabun*100, na.rm=T))  %>%
    mutate(total=str_c(round(total_mean,2),"±",round(total_sd,2))) %>% 
    arrange(-total_mean) %>% 
    dplyr::select(phylum,total, Population) %>% 
    tt()

phylum	total	Population
p__Bacteroidota	41.93±15.4	Cold_wet
p__Bacillota_A	33.96±16.7	Cold_wet
p__Bacillota	5.73±4.73	Cold_wet
p__Campylobacterota	5.39±4.64	Cold_wet
p__Verrucomicrobiota	4.36±3.06	Cold_wet
p__Pseudomonadota	3.75±3.76	Cold_wet
p__Desulfobacterota	2.31±1.6	Cold_wet
p__Bacillota_C	0.85±0.82	Cold_wet
p__Cyanobacteriota	0.64±0.83	Cold_wet
p__Bacillota_B	0.5±0.77	Cold_wet
p__Fusobacteriota	0.37±1.28	Cold_wet
p__Spirochaetota	0.12±0.35	Cold_wet
p__Actinomycetota	0.09±0.22	Cold_wet
p__Chlamydiota	0±0	Cold_wet
p__Elusimicrobiota	0±0	Cold_wet

#Hot and dry population
phylum_summary %>%
    left_join(sample_metadata, by = join_by(sample == Tube_code))  %>%
    group_by(phylum, Population) %>%
    filter(Population=="Hot_dry")%>%
    filter(time_point=="0_Wild") %>%
    summarise(total_mean=mean(relabun*100, na.rm=T),
              total_sd=sd(relabun*100, na.rm=T))  %>%
    mutate(total=str_c(round(total_mean,2),"±",round(total_sd,2))) %>% 
    arrange(-total_mean) %>% 
    dplyr::select(phylum,total, Population) %>% 
    tt()

phylum	total	Population
p__Bacillota_A	33.9±16.17	Hot_dry
p__Bacteroidota	27.63±17.64	Hot_dry
p__Campylobacterota	13.53±20.98	Hot_dry
p__Pseudomonadota	10.92±9.94	Hot_dry
p__Bacillota	6.09±7.66	Hot_dry
p__Spirochaetota	2.66±2.78	Hot_dry
p__Desulfobacterota	1.69±1.85	Hot_dry
p__Fusobacteriota	1.34±1.73	Hot_dry
p__Bacillota_B	0.76±0.73	Hot_dry
p__Bacillota_C	0.73±0.5	Hot_dry
p__Cyanobacteriota	0.52±0.65	Hot_dry
p__Chlamydiota	0.12±0.18	Hot_dry
p__Verrucomicrobiota	0.06±0.1	Hot_dry
p__Elusimicrobiota	0.05±0.14	Hot_dry
p__Actinomycetota	0±0	Hot_dry

phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal() + 
        theme(legend.position="none") +
        labs(y="Phylum",x="Relative abundance")

5.2 Taxonomy boxplot

5.2.1 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()

family	mean	sd
f__Bacteroidaceae	2.260146e-01	0.1384706235
f__Lachnospiraceae	1.410833e-01	0.1062953893
f__Tannerellaceae	1.045659e-01	0.0799894745
f__Helicobacteraceae	5.448546e-02	0.0937279764
f__Mycoplasmoidaceae	3.756572e-02	0.0767893776
f__Erysipelotrichaceae	3.536287e-02	0.0452140267
f__UBA3700	3.456595e-02	0.0565495010
f__Marinifilaceae	2.794365e-02	0.0272474083
f__Rikenellaceae	2.725202e-02	0.0471735191
f__Enterobacteriaceae	2.687327e-02	0.0929254305
f__Coprobacillaceae	2.627823e-02	0.0907456387
f__	2.465083e-02	0.0781013121
f__Desulfovibrionaceae	2.373771e-02	0.0367718116
f__DTU072	2.183007e-02	0.0377159876
f__Ruminococcaceae	1.832093e-02	0.0428115194
f__Rhizobiaceae	1.579679e-02	0.0779688169
f__LL51	1.556592e-02	0.0616955422
f__UBA3830	1.512118e-02	0.0441855701
f__Akkermansiaceae	1.224165e-02	0.0317653885
f__Chlamydiaceae	1.096188e-02	0.0607502761
f__Fusobacteriaceae	1.055726e-02	0.0286383947
f__CAG-239	9.138683e-03	0.0152490206
f__Enterococcaceae	8.437601e-03	0.0473906561
f__Gastranaerophilaceae	7.848357e-03	0.0146292952
f__Oscillospiraceae	6.624721e-03	0.0075288565
f__UBA1997	6.613196e-03	0.0315103296
f__Streptococcaceae	6.600789e-03	0.0348230465
f__UBA1242	4.266475e-03	0.0147768596
f__Brevinemataceae	4.098862e-03	0.0125062564
f__Acutalibacteraceae	3.498766e-03	0.0111374416
f__RUG11792	2.921450e-03	0.0255676374
f__Clostridiaceae	2.855351e-03	0.0174153876
f__UBA660	2.600647e-03	0.0118148140
f__Peptococcaceae	2.556218e-03	0.0049947786
f__Acidaminococcaceae	1.980431e-03	0.0051045211
f__CAG-508	1.874529e-03	0.0065256902
f__MGBC116941	1.783700e-03	0.0077801927
f__Moraxellaceae	1.540093e-03	0.0099192011
f__RUG14156	1.428152e-03	0.0045670616
f__Staphylococcaceae	1.411833e-03	0.0051727635
f__Anaerovoracaceae	1.410739e-03	0.0027876311
f__Elusimicrobiaceae	1.088209e-03	0.0062780187
f__CAG-288	9.840222e-04	0.0061275213
f__Anaerotignaceae	9.320656e-04	0.0041174457
f__CALVMC01	7.793540e-04	0.0044385554
f__Eggerthellaceae	6.643755e-04	0.0021620275
f__Massilibacillaceae	6.322621e-04	0.0016561037
f__Mycobacteriaceae	6.168531e-04	0.0061497354
f__UBA1820	4.705627e-04	0.0013078764
f__CAG-274	4.686117e-04	0.0022415212
f__Arcobacteraceae	3.928587e-04	0.0050156837
f__Burkholderiaceae_C	3.835606e-04	0.0048969735
f__Muribaculaceae	3.508548e-04	0.0009525792
f__UBA932	3.295199e-04	0.0011408058
f__Hepatoplasmataceae	3.099135e-04	0.0039567109
f__Rhodobacteraceae	3.068016e-04	0.0039169801
f__Weeksellaceae	2.873650e-04	0.0032049404
f__Eubacteriaceae	1.707442e-04	0.0006844943
f__Sphingobacteriaceae	1.561202e-04	0.0012685229
f__Devosiaceae	1.544841e-04	0.0015368528
f__Pumilibacteraceae	1.324439e-04	0.0007783049
f__WRAU01	9.956857e-05	0.0012712064
f__Peptostreptococcaceae	2.371535e-05	0.0003027773

family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

5.2.2 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()

genus	mean	sd
g__Bacteroides	1.374441e-01	0.0923562611
g__Parabacteroides	9.843371e-02	0.0803813676
g__Phocaeicola	7.109518e-02	0.0799972428
g__Helicobacter_J	3.115738e-02	0.0603033642
g__Mycoplasmoides	3.115302e-02	0.0765633496
g__Odoribacter	2.605667e-02	0.0268723184
g__Roseburia	2.425344e-02	0.0567044550
g__NHYM01	2.332808e-02	0.0810376677
g__Alistipes	2.224698e-02	0.0287419149
g__Coprobacillus	2.070109e-02	0.0894282233
g__Agrobacterium	1.579679e-02	0.0779688169
g__Akkermansia	1.224165e-02	0.0317653885
g__Fusobacterium_A	1.046073e-02	0.0286441117
g__Kineothrix	9.107908e-03	0.0416218616
g__Proteus	8.976683e-03	0.0694135711
g__Dielma	8.687357e-03	0.0090713197
g__CAG-95	8.238073e-03	0.0207930753
g__JAAYNV01	7.265789e-03	0.0179169564
g__Desulfovibrio	7.219276e-03	0.0214990147
g__UBA866	7.016767e-03	0.0295145125
g__Enterococcus	6.966137e-03	0.0463712943
g__Lactococcus	6.600789e-03	0.0348230465
g__Ureaplasma	6.412700e-03	0.0139267552
g__Parabacteroides_B	6.132159e-03	0.0101543965
g__Lacrimispora	6.028411e-03	0.0098068179
g__CALXRO01	5.977964e-03	0.0313982647
g__Citrobacter	5.896711e-03	0.0340533686
g__NSJ-61	5.745781e-03	0.0202234473
g__Breznakia	5.530147e-03	0.0240721461
g__Clostridium_AQ	5.522246e-03	0.0123487838
g__Bilophila	5.044501e-03	0.0089558435
g__Hungatella_A	4.964136e-03	0.0096921078
g__Escherichia	4.342538e-03	0.0270859242
g__Salmonella	4.319018e-03	0.0148769561
g__UMGS1251	4.312965e-03	0.0073601071
g__MGBC136627	4.305492e-03	0.0164533523
g__Hungatella	4.150386e-03	0.0194068227
g__Clostridium_Q	4.146767e-03	0.0052575243
g__Brevinema	4.098862e-03	0.0125062564
g__Thomasclavelia	4.046233e-03	0.0110779301
g__Scatousia	3.752075e-03	0.0104403539
g__Mailhella	3.745039e-03	0.0104110785
g__Copromonas	3.643508e-03	0.0050495456
g__Enterocloster	3.613702e-03	0.0047492729
g__Ventrimonas	3.566172e-03	0.0071931788
g__Fournierella	3.313097e-03	0.0063192740
g__Limenecus	3.230504e-03	0.0066725343
g__Mucinivorans	3.006847e-03	0.0379999623
g__Lawsonia	2.916613e-03	0.0103686789
g__MGBC133411	2.902785e-03	0.0074333461
g__Caccovivens	2.887473e-03	0.0112659902
g__Sarcina	2.855351e-03	0.0174153876
g__Eisenbergiella	2.796704e-03	0.0069384489
g__Bacteroides_G	2.781473e-03	0.0352463088
g__CAJLXD01	2.730769e-03	0.0088951735
g__Acetatifactor	2.654208e-03	0.0055286194
g__Blautia	2.598789e-03	0.0062369300
g__Velocimicrobium	2.235984e-03	0.0067748392
g__C-19	2.235603e-03	0.0048296119
g__CAZU01	2.189719e-03	0.0066369837
g__Negativibacillus	2.145239e-03	0.0056002700
g__Intestinimonas	2.003816e-03	0.0035552824
g__Rikenella	1.998193e-03	0.0037323264
g__Phascolarctobacterium	1.980431e-03	0.0051045211
g__Butyricimonas	1.886974e-03	0.0042483569
g__RGIG6463	1.855727e-03	0.0040258495
g__MGBC116941	1.783700e-03	0.0077801927
g__JALFVM01	1.712574e-03	0.0038669765
g__Oscillibacter	1.546231e-03	0.0025273862
g__Acinetobacter	1.540093e-03	0.0099192011
g__Pseudoflavonifractor	1.489286e-03	0.0027026675
g__Citrobacter_A	1.444197e-03	0.0061395639
g__Staphylococcus	1.411833e-03	0.0051727635
g__14-2	1.228249e-03	0.0098038984
g__RGIG4733	1.225656e-03	0.0038271024
g__Beduini	1.217163e-03	0.0025500013
g__Scatocola	1.162463e-03	0.0045748144
g__Enterococcus_A	1.123893e-03	0.0100947603
g__UBA1436	1.088209e-03	0.0062780187
g__Faecisoma	1.057958e-03	0.0056043491
g__RGIG9287	9.993630e-04	0.0094920364
g__CAG-345	9.840222e-04	0.0061275213
g__Lachnotalea	9.593025e-04	0.0033990676
g__Blautia_A	9.546974e-04	0.0029560949
g__Ruthenibacterium	8.602962e-04	0.0024818365
g__CAG-269	8.276280e-04	0.0048017072
g__Marseille-P3106	8.230475e-04	0.0017874580
g__WRHT01	6.666234e-04	0.0027445999
g__Eggerthella	6.643755e-04	0.0021620275
g__CHH4-2	6.371240e-04	0.0020328940
g__Corynebacterium	6.168531e-04	0.0061497354
g__Serratia_A	6.076344e-04	0.0077577570
g__Anaerotruncus	6.058602e-04	0.0016447558
g__RUG14156	5.735678e-04	0.0021869659
g__RGIG1896	5.683407e-04	0.0051791669
g__IOR16	5.574841e-04	0.0016418264
g__Faecimonas	5.146607e-04	0.0054508437
g__CAG-56	5.096368e-04	0.0016613952
g__MGBC140009	4.851579e-04	0.0024491799
g__CALURL01	4.805911e-04	0.0017020401
g__Merdimorpha	4.705627e-04	0.0013078764
g__RGIG8482	4.560993e-04	0.0030287706
g__Enterobacter	4.223379e-04	0.0042068345
g__Klebsiella	4.203682e-04	0.0049802041
g__Caccenecus	4.086273e-04	0.0018112589
g__Aliarcobacter	3.928587e-04	0.0050156837
g__Scatenecus	3.851876e-04	0.0018282510
g__Alcaligenes	3.835606e-04	0.0048969735
g__Plesiomonas	3.766988e-04	0.0027593254
g__JAHHSE01	3.529590e-04	0.0014998851
g__HGM05232	3.508548e-04	0.0009525792
g__Enterococcus_B	3.475714e-04	0.0022665993
g__Egerieousia	3.295199e-04	0.0011408058
g__Stoquefichus	3.137462e-04	0.0020871798
g__Hepatoplasma	3.099135e-04	0.0039567109
g__Paracoccus	3.068016e-04	0.0039169801
g__Moheibacter	2.873650e-04	0.0032049404
g__Scatomorpha	2.738230e-04	0.0010358302
g__Emergencia	2.601331e-04	0.0013298673
g__UBA7185	2.523935e-04	0.0014817660
g__Eubacterium	1.707442e-04	0.0006844943
g__Sphingobacterium	1.561202e-04	0.0012685229
g__Devosia	1.544841e-04	0.0015368528
g__Anaerosporobacter	1.507638e-04	0.0012978048
g__Caccomorpha	1.434035e-04	0.0010730603
g__UBA2658	1.355578e-04	0.0007332702
g__Protoclostridium	1.324439e-04	0.0007783049
g__Angelakisella	1.315171e-04	0.0009387427
g__Cetobacterium_A	9.652924e-05	0.0008876688
g__Rahnella	6.708891e-05	0.0008565338
g__Peptostreptococcus	2.371535e-05	0.0003027773

genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

6.3.1.2 Number of samples used

[1] 27

beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

6.3.1.3 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.975
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.979
Podarcis_muralis            0.97302                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.542719	0.2095041	6.625717	0.001
Residual	25	5.820951	0.7904959	NA	NA
Total	26	7.363669	1.0000000	NA	NA

6.3.1.4 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999  0.948
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.946
Podarcis_muralis            0.94763                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.918266	0.2608511	8.822682	0.001
Residual	25	5.435610	0.7391489	NA	NA
Total	26	7.353876	1.0000000	NA	NA

6.3.1.5 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.135
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.133
Podarcis_muralis            0.12705                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.3218613	0.2162815	6.899207	0.001
Residual	25	1.1662981	0.7837185	NA	NA
Total	26	1.4881594	1.0000000	NA	NA

6.3.1.6 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.018367 0.018367 1.5597    999  0.203
Residuals 25 0.294402 0.011776                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.209
Podarcis_muralis            0.22328                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0858578	0.172879	5.225323	0.053
Residual	25	0.4107775	0.827121	NA	NA
Total	26	0.4966352	1.000000	NA	NA

beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

accli <- meta %>%
  filter(time_point == "1_Acclimation")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))

accli_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation")

6.3.2.1 Number of samples used

[1] 27

beta_div_richness_accli<-hillpair(data=accli.counts, q=0)
beta_div_neutral_accli<-hillpair(data=accli.counts, q=1)
beta_div_phylo_accli<-hillpair(data=accli.counts, q=1, tree=genome_tree)
beta_div_func_accli<-hillpair(data=accli.counts, q=1, dist=dist)

6.3.2.2 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.11796 0.117959 12.963    999  0.003 **
Residuals 25 0.22748 0.009099                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.002
Hot_dry  0.0013711        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.639807	0.179834	5.481634	0.001
Residual	25	7.478640	0.820166	NA	NA
Total	26	9.118447	1.000000	NA	NA

6.3.2.3 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07844 0.078443 5.2384    999  0.032 *
Residuals 25 0.37437 0.014975                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.029
Hot_dry  0.030815        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.947003	0.2306127	7.493387	0.001
Residual	25	6.495736	0.7693873	NA	NA
Total	26	8.442739	1.0000000	NA	NA

6.3.2.4 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.06739 0.067395 2.9532    999  0.104
Residuals 25 0.57052 0.022821                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.095
Hot_dry  0.098068        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.2441653	0.1224638	3.488854	0.032
Residual	25	1.7496100	0.8775362	NA	NA
Total	26	1.9937754	1.0000000	NA	NA

6.3.2.5 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02496 0.024955 0.6729    999  0.454
Residuals 25 0.92714 0.037085                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.448
Hot_dry   0.41979        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.0279454	0.0248037	0.6358634	0.466
Residual	25	1.0987171	0.9751963	NA	NA
Total	26	1.1266624	1.0000000	NA	NA

beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

accli1 <- meta  %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))

accli1_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

6.3.3.1 Number of samples used

[1] 54

beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)

6.3.3.1.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05014 0.050145 6.2252    999  0.011 *
Residuals 52 0.41886 0.008055                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                0_Wild 1_Acclimation
0_Wild                         0.014
1_Acclimation 0.015808              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.6172653	0.0360987	3.933397	0.001
species	1	2.8279677	0.1653842	18.020647	0.001
individual	25	9.5739861	0.5599025	2.440331	0.001
Residual	26	4.0801621	0.2386146	NA	NA
Total	53	17.0993812	1.0000000	NA	NA

6.3.3.1.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0199 0.0199035 2.1213    999  0.139
Residuals 52 0.4879 0.0093827                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.142
1_Acclimation 0.15128              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.9050519	0.0541893	6.651487	0.001
species	1	3.3236300	0.1989999	24.426315	0.001
individual	25	8.9352276	0.5349902	2.626702	0.001
Residual	26	3.5377576	0.2118206	NA	NA
Total	53	16.7016671	1.0000000	NA	NA

6.3.3.1.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01334 0.013340 0.6524    999  0.435
Residuals 52 1.06332 0.020449                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.437
1_Acclimation 0.42294              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.2890434	0.0766494	7.532050	0.001
species	1	0.3508889	0.0930498	9.143655	0.001
individual	25	2.1332925	0.5657133	2.223620	0.001
Residual	26	0.9977533	0.2645874	NA	NA
Total	53	3.7709782	1.0000000	NA	NA

6.3.3.1.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0123 0.012300 0.4817    999  0.482
Residuals 52 1.3277 0.025533                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.487
1_Acclimation 0.49073              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.0448774	0.0269021	2.355512	0.163
species	1	0.0973005	0.0583275	5.107077	0.034
individual	25	1.0306426	0.6178264	2.163841	0.055
Residual	26	0.4953546	0.2969440	NA	NA
Total	53	1.6681751	1.0000000	NA	NA

beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

6.3.4 4. Do the antibiotics work?

6.3.4.1 Acclimation vs antibiotics

treat <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

6.3.4.2 Number of samples used

[1] 50

beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)

6.3.4.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.025318 0.0253178 6.021    999  0.023 *
Residuals 48 0.201837 0.0042049                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.027
2_Antibiotics      0.017817              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.888584	0.0949462	6.361098	0.001
species	1	2.117109	0.1064350	7.130814	0.001
individual	25	9.353701	0.4702455	1.260199	0.005
Residual	22	6.531709	0.3283734	NA	NA
Total	49	19.891103	1.0000000	NA	NA

6.3.4.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.039587 0.039587 6.8387    999  0.013 *
Residuals 48 0.277854 0.005789                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.015
2_Antibiotics      0.011886              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	2.024181	0.1063620	9.051981	0.001
species	1	2.853103	0.1499183	12.758858	0.001
individual	25	9.234189	0.4852168	1.651783	0.001
Residual	22	4.919584	0.2585029	NA	NA
Total	49	19.031057	1.0000000	NA	NA

6.3.4.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.58372 0.58372 35.413    999  0.001 ***
Residuals 48 0.79119 0.01648                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.001
2_Antibiotics    2.9795e-07              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.8065206	0.2113909	18.636551	0.001
species	1	0.7903334	0.0924813	8.153292	0.001
individual	25	3.8164689	0.4465860	1.574869	0.005
Residual	22	2.1325541	0.2495419	NA	NA
Total	49	8.5458771	1.0000000	NA	NA

6.3.4.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.18591 0.185914 5.0679    999  0.023 *
Residuals 48 1.76088 0.036685                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.032
2_Antibiotics      0.028989              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.8020952	0.3750193	33.6195614	0.001
species	1	0.0031247	0.0006503	0.0582938	0.848
individual	25	1.8208629	0.3789249	1.3587875	0.222
Residual	22	1.1792568	0.2454055	NA	NA
Total	49	4.8053396	1.0000000	NA	NA

beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

6.3.5 5. Does the FMT work?

6.3.5.1 Comparison between FMT2 vs Post-FMT2

#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

transplant3<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")%>%
  column_to_rownames("newID")

transplant3_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

full_counts<-temp_genome_counts %>%
    t()%>%
    as.data.frame()%>%
    rownames_to_column("Tube_code")%>%
    full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))

transplant3_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

identical(sort(colnames(transplant3_counts)),sort(as.character(rownames(transplant3))))

6.3.5.2 Number of samples used

[1] 49

beta_div_richness_transplant3<-hillpair(data=transplant3_counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3_counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3_counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3_counts, q=1, dist=dist)

#Arrange of metadata dataframe
transplant3_arrange<-transplant3[labels(beta_div_neutral_transplant3$S),]

6.3.5.2.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.180473	0.0855095	6.984555	0.001
time_point	1	0.860906	0.0623612	5.093759	0.001
type	1	1.459433	0.1057165	8.635089	0.001
individual	24	6.755100	0.4893170	1.665341	0.001
Residual	21	3.549250	0.2570959	NA	NA
Total	48	13.805162	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.4169018 5.739828 0.15622903   0.001      0.003   *
2   Control vs Hot_control  1 2.0940966 8.509112 0.21005427   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3004618 1.265034 0.04179854   0.159      0.477    

6.3.5.2.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.2800927	0.0939787	8.796453	0.001
time_point	1	0.9350566	0.0686477	6.425458	0.001
type	1	1.9135997	0.1404879	13.149743	0.001
individual	24	6.4363516	0.4725281	1.842870	0.001
Residual	21	3.0559984	0.2243577	NA	NA
Total	48	13.6210990	1.0000000	NA	NA

                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.8758788  8.282671 0.21084796   0.001      0.003   *
2   Control vs Hot_control  1 2.4396317 10.635546 0.24945256   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3158428  1.394345 0.04587515   0.126      0.378    

6.3.5.2.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.1400466	0.0952654	6.956436	0.002
time_point	1	0.1138047	0.0774145	5.652939	0.001
type	1	0.1432667	0.0974558	7.116383	0.001
individual	24	0.6501795	0.4422784	1.345663	0.041
Residual	21	0.4227709	0.2875859	NA	NA
Total	48	1.4700683	1.0000000	NA	NA

                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.14387705 5.735321 0.15612552   0.001      0.003   *
2   Control vs Hot_control  1 0.22715701 9.044894 0.22036587   0.001      0.003   *
3 Treatment vs Hot_control  1 0.04648319 1.704277 0.05550617   0.129      0.387    

6.3.5.2.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0092808	0.0077189	0.4182529	0.493
time_point	1	-0.0061674	-0.0051295	-0.2779456	0.895
type	1	0.0831052	0.0691191	3.7452726	0.093
individual	24	0.6501528	0.5407359	1.2208414	0.345
Residual	21	0.4659767	0.3875556	NA	NA
Total	48	1.2023481	1.0000000	NA	NA

                     pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1     Control vs Treatment  1  0.078539743  4.59293783  0.129040706   0.069      0.207    
2   Control vs Hot_control  1  0.052468954  2.13675422  0.062593948   0.165      0.495    
3 Treatment vs Hot_control  1 -0.002340352 -0.07432315 -0.002569452   0.887      1.000    

beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

6.3.5.3 Comparison between the different experimental time points

6.3.5.3.1 Comparison against Wild samples

The estimated time for calculating the 2850 pairwise combinations is 14 seconds.

ggarrange( p1, p2, p3, ncol=3, nrow=1, common.legend = TRUE, legend="right")

6.3.5.3.2 Comparison against Acclimation samples

The estimated time for calculating the 2850 pairwise combinations is 12 seconds.

ggarrange( p1, p2, p3, ncol=3, nrow=1, common.legend = TRUE, legend="right")

6.3.5.3.3 Comparison against Acclimation samples with combined Transplant samples

The estimated time for calculating the 5151 pairwise combinations is 22 seconds.

ggarrange( p1, p2, p3, ncol=3, nrow=1, common.legend = TRUE, legend="right")

6.3.6 6. Are there differences between the control and the treatment group?

6.3.6.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

6.3.6.2 Number of samples used

[1] 26

beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]

6.3.6.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.017675 0.0088373 2.3825    999  0.099 .
Residuals 23 0.085312 0.0037092                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0060000     0.661
Hot_control 0.0068795                 0.212
Treatment   0.6248469   0.2084296          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.6340254	0.0768024	2.065607	0.004
type	1	0.5615418	0.0680222	1.829461	0.010
Residual	23	7.0597099	0.8551754	NA	NA
Total	25	8.2552771	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.5615418 1.729004 0.1033537   0.016      0.048   .
2   Control vs Hot_control  1 0.8438429 2.793772 0.1486541   0.002      0.006   *
3 Treatment vs Hot_control  1 0.3734921 1.268929 0.0779971   0.098      0.294    

6.3.6.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999  0.564
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.21500     0.957
Hot_control 0.21166                 0.467
Treatment   0.95468     0.43604          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.7907904	0.1076445	3.056657	0.001
type	1	0.6051778	0.0823784	2.339205	0.009
Residual	23	5.9503501	0.8099772	NA	NA
Total	25	7.3463184	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.6051778 2.250849 0.13047758   0.009      0.027   .
2   Control vs Hot_control  1 1.0528902 4.143637 0.20570451   0.001      0.003   *
3 Treatment vs Hot_control  1 0.4150076 1.637268 0.09840968   0.046      0.138    

6.3.6.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999   0.91
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92400     0.694
Hot_control 0.91505                 0.803
Treatment   0.63312     0.73046          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0560850	0.0531376	1.3149967	0.280
type	1	0.0184254	0.0174571	0.4320099	0.791
Residual	23	0.9809570	0.9294053	NA	NA
Total	25	1.0554673	1.0000000	NA	NA

                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01842535 0.4144162 0.02688498   0.775      1.000    
2   Control vs Hot_control  1 0.05987967 1.7387847 0.09802164   0.117      0.351    
3 Treatment vs Hot_control  1 0.03212966 0.6477782 0.04139746   0.692      1.000    

6.3.6.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.00400 0.0020014 0.145    999  0.864
Residuals 23 0.31753 0.0138057                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.59900     0.750
Hot_control 0.59817                 0.849
Treatment   0.75141     0.83718          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0024979	0.0033024	0.0900845	0.638
type	1	0.1161466	0.1535542	4.1887855	0.063
Residual	23	0.6377435	0.8431434	NA	NA
Total	25	0.7563879	1.0000000	NA	NA

                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.11614656 4.724791 0.23953568   0.069      0.207    
2   Control vs Hot_control  1 0.05000930 1.704826 0.09629160   0.224      0.672    
3 Treatment vs Hot_control  1 0.01235859 0.423812 0.02747777   0.503      1.000    

beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.6.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

6.3.6.4 Number of samples used

[1] 27

beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]

6.3.6.4.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999   0.83
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.70300     0.786
Hot_control 0.67789                 0.624
Treatment   0.79246     0.59820          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	1.504341	0.1967776	2.939822	0.001
Residual	24	6.140538	0.8032224	NA	NA
Total	26	7.644879	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 0.6463814 2.560441 0.1379515   0.001      0.003   *
2 Treatment vs Hot_control  1 0.4796256 1.916520 0.1069694   0.001      0.003   *
3   Control vs Hot_control  1 1.1305044 4.268317 0.2105906   0.001      0.003   *

6.3.6.4.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.468
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.46700     0.694
Hot_control 0.44675                 0.241
Treatment   0.65989     0.25095          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	1.923807	0.2603795	4.224537	0.001
Residual	24	5.464666	0.7396205	NA	NA
Total	26	7.388473	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 1.0227481 4.648335 0.2251191   0.001      0.003   *
2 Treatment vs Hot_control  1 0.5010202 2.206532 0.1211945   0.001      0.003   *
3   Control vs Hot_control  1 1.3619424 5.771031 0.2650785   0.001      0.003   *

6.3.6.4.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.959
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92800     0.874
Hot_control 0.93765                 0.769
Treatment   0.83933     0.76015          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	0.1594363	0.2042241	3.079623	0.001
Residual	24	0.6212564	0.7957759	NA	NA
Total	26	0.7806927	1.0000000	NA	NA

                     pairs Df  SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 0.05927454 2.382025 0.1295845   0.021      0.063    
2 Treatment vs Hot_control  1 0.06906280 2.722460 0.1454115   0.005      0.015   .
3   Control vs Hot_control  1 0.11081709 4.043656 0.2017424   0.001      0.003   *

6.3.6.4.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01126 0.0056302 0.2861    999  0.773
Residuals 24 0.47233 0.0196806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.54800     0.634
Hot_control 0.48255                 0.790
Treatment   0.60116     0.75643          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	-0.0038724	-0.0056213	-0.0670788	0.936
Residual	24	0.6927468	1.0056213	NA	NA
Total	26	0.6888744	1.0000000	NA	NA

                     pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1     Treatment vs Control  1 -0.008527330 -0.46290555 -0.029793572   0.853          1    
2 Treatment vs Hot_control  1 -0.001648721 -0.04717131 -0.002956924   0.910          1    
3   Control vs Hot_control  1  0.004367477  0.13147026  0.008149924   0.690          1    

beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

7 Functional differences

load("data/data.Rdata")

7.1 Data preparation

# Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts, GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts$genome, ]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
  select_if(~ !is.numeric(.) || sum(.) != 0)

elements <- GIFTs_elements_filtered %>%
  as.data.frame()

# Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered, GIFT_db)
functions <- GIFTs_functions %>%
  as.data.frame()

# Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions, GIFT_db)
domains <- GIFTs_domains %>%
  as.data.frame()

# Get community-weighed average GIFTs per sample
GIFTs_elements_community <- to.community(GIFTs_elements_filtered, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db)

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

7.2 Genomes GIFT profiles

GIFTs_elements %>%
  as_tibble(., rownames = "MAG") %>%
  reshape2::melt() %>%
  rename(Code_element = variable, GIFT = value) %>%
  inner_join(GIFT_db,by="Code_element") %>%
  ggplot(., aes(x=Code_element, y=MAG, fill=GIFT, group=Code_function))+
    geom_tile()+
    scale_y_discrete(guide = guide_axis(check.overlap = TRUE))+
    scale_x_discrete(guide = guide_axis(check.overlap = TRUE))+
    scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
    facet_grid(. ~ Code_function, scales = "free", space = "free")+
    theme_grey(base_size=8)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text.x = element_text(angle = 90))

7.3 Function level

GIFTs_functions_community %>%
    as.data.frame() %>%
    rownames_to_column(var="sample") %>%
    filter(sample!="AD69") %>%
    pivot_longer(!sample,names_to="trait",values_to="gift") %>%
    left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
    ggplot(aes(x=trait,y=time_point,fill=gift)) +
        geom_tile(colour="white", size=0.2)+
        scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
        facet_grid(type ~ ., scales="free",space="free")

7.4 Element level

GIFTs_elements_community_merged<-GIFTs_elements_community %>%
    as.data.frame() %>%
    rownames_to_column(var="sample") %>%
    filter(sample!="AD69") %>%
    pivot_longer(!sample,names_to="trait",values_to="gift") %>%
    left_join(sample_metadata, by = join_by(sample == Tube_code))%>%
    mutate(functionid = substr(trait, 1, 3)) %>%
    mutate(trait = case_when(
      trait %in% GIFT_db$Code_element ~ GIFT_db$Element[match(trait, GIFT_db$Code_element)],
      TRUE ~ trait
    )) %>%
    mutate(functionid = case_when(
      functionid %in% GIFT_db$Code_function ~ GIFT_db$Function[match(functionid, GIFT_db$Code_function)],
      TRUE ~ functionid
    )) %>%
    mutate(trait=factor(trait,levels=unique(GIFT_db$Element))) %>%
    mutate(functionid=factor(functionid,levels=unique(GIFT_db$Function)))

# Create an interaction variable for time_point and sample
GIFTs_elements_community_merged$interaction_var <- interaction(GIFTs_elements_community_merged$sample, GIFTs_elements_community_merged$time_point)
  
ggplot(GIFTs_elements_community_merged,aes(x=interaction_var,y=trait,fill=gift)) +
        geom_tile(colour="white", linewidth=0.2)+
        scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
        facet_grid(functionid ~ type, scales="free",space="free") +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=5),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Time_point",fill="GIFT")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x))

7.5 Comparison of samples from the 0 Time_point (0_Wild)

7.5.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.327 0.0244
2 Podarcis_muralis  0.346 0.0194

7.5.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  select(c(1:21, 27)) %>%
  pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
  ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
    geom_boxplot() +
    scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c","#76b183")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.5.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.375 0.0315
2 Podarcis_muralis  0.390 0.0208

7.5.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.313 0.0329
2 Podarcis_muralis  0.345 0.0233

sample_metadata_wild <- sample_metadata%>% 
  filter(time_point == "0_Wild")

element_gift_wild <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_wild[c(1,3)], by="Tube_code")

# Find numeric columns
numeric_cols <- sapply(element_gift_wild, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_wild[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_wild[, !numeric_cols | colnames(element_gift_wild) %in% nonzero_numeric_cols]

significant_elements_wild <- filtered_data %>%
  pivot_longer(-c(Tube_code,species), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ species, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_wild  %>% 
  dplyr::select(-c(species))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_wild$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_wild[c(1,3)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(species)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-species)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_wild[c(1,4)], by = join_by(Tube_code == Tube_code))

colNames <- names(element_gift_names)[2:30] #always check names(element_gift_names) first to know where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=Population, y=.data[[i]], color = Population)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.6 Comparison of samples from the 1st Time_point (1_Acclimation)

7.6.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.348 0.0158
2 Podarcis_muralis  0.331 0.0321

7.6.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  select(c(1:21, 27)) %>%
  pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
  ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
    geom_boxplot() +
    scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c","#76b183")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.6.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.395 0.0211
2 Podarcis_muralis  0.370 0.0307

7.6.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.350 0.0225
2 Podarcis_muralis  0.332 0.0316

sample_metadata_accli <- sample_metadata%>% 
  filter(time_point == "1_Acclimation")

element_gift_accli <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_accli[c(1,3)], by="Tube_code")

# Find numeric columns
numeric_cols <- sapply(element_gift_accli, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_accli[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_accli[, !numeric_cols | colnames(element_gift_accli) %in% nonzero_numeric_cols]

significant_elements_accli <- filtered_data %>%
  pivot_longer(-c(Tube_code,species), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ species, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_accli  %>% 
  dplyr::select(-c(species))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_accli$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_accli[c(1,3)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(species)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-species)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_accli[c(1,4)], by = join_by(Tube_code == Tube_code))

colNames <- names(element_gift_names)[2:14] #always check names(element_gift_names) first to know where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=Population, y=.data[[i]], color = Population)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.7 Comparison of samples from the 5th Time_point (5_Post-FMT1)

7.7.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.373 0.0247
2 Hot_control 0.372 0.0367
3 Treatment   0.353 0.0186

7.7.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  select(c(1:21, 27)) %>%
  pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
  ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
    geom_boxplot() +
    scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c","#76b183")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.7.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.412 0.0213
2 Hot_control 0.415 0.0437
3 Treatment   0.391 0.0263

7.7.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.380 0.0280
2 Hot_control 0.379 0.0372
3 Treatment   0.359 0.0214

sample_metadata_tm5 <- sample_metadata%>% 
  filter(time_point == "5_Post-FMT1")%>% 
  filter(type != "Hot_control")

element_gift_tm5 <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(sample_metadata_tm5 %>% select(1, 7), by = "Tube_code")

# Find numeric columns
numeric_cols <- sapply(element_gift_tm5, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_tm5[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_tm5[, !numeric_cols | colnames(element_gift_tm5) %in% nonzero_numeric_cols]

significant_elements_tm5 <- filtered_data %>%
  pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ type, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_value < 0.05)  %>% #take into account that p_value is used and not p_adjust
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_tm5  %>% 
  dplyr::select(-c(type))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_tm5$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_tm5[c(1,7)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(type)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-type)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_tm5[c(1,7)], by = join_by(Tube_code == Tube_code))

colNames <- names(element_gift_names)[2:14] #always check names(element_gift_names) first to now where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=type, y=.data[[i]], color = type)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.8 Comparison of samples from the 6th Time_point (6_Post-FMT2)

7.8.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.352 0.0223
2 Hot_control 0.350 0.0293
3 Treatment   0.346 0.0255

7.8.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  select(c(1:21, 27)) %>%
  pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db$Code_function ~ GIFT_db$Function[match(trait, GIFT_db$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db$Function))) %>%
  ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
    geom_boxplot() +
    scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c","#76b183")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.8.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.399 0.0171
2 Hot_control 0.388 0.0321
3 Treatment   0.392 0.0240

7.8.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.357 0.0215
2 Hot_control 0.347 0.0302
3 Treatment   0.350 0.0293

sample_metadata_TM6 <- sample_metadata%>% 
  filter(time_point == "6_Post-FMT2")%>% 
  filter(type != "Hot_control")

element_gift_TM6 <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(sample_metadata_TM6 %>% select(1, 7), by = "Tube_code")

# Find numeric columns
numeric_cols <- sapply(element_gift_TM6, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_TM6[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_TM6[, !numeric_cols | colnames(element_gift_TM6) %in% nonzero_numeric_cols]

significant_elements_TM6 <- filtered_data %>%
  pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ type, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_value < 0.05)  %>% #take into account that p_value is used and not p_adjust
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_TM6  %>% 
  dplyr::select(-c(type))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_TM6$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_TM6[c(1,7)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(type)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-type)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_TM6[c(1,7)], by = join_by(Tube_code == Tube_code))

colNames <- names(element_gift_names)[2:20] #always check names(element_gift_names) first to now where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=type, y=.data[[i]], color = type)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.9 Domain level

7.9.1 Comparison of samples from the 0 Time_point (0_Wild)

#Merge the functional domains with the metadata
merge_gift_wild<- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_wild, by="Tube_code")

#Biosynthesis
p1 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Degradation
p2 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Structure
p3 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Structure,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Overall
p4 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Overall,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

7.9.2 Comparison of samples from the 1st Time_point (1_Acclimation)

#Merge the functional domains with the metadata
merge_gift_accli<- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_accli, by="Tube_code")

#Biosynthesis
p1 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Degradation
p2 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Structure
p3 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Structure,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Overall
p4 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Overall,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

7.9.3 Comparison of samples from the 5th Time_point (5_Post-FMT1)

#Merge the functional domains with the metadata
merge_gift_tm5<- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_tm5, by="Tube_code")

#Biosynthesis
p1 <-merge_gift_tm5 %>%
  ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Degradation
p2 <-merge_gift_tm5 %>%
  ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Structure
p3 <-merge_gift_tm5 %>%
  ggplot(aes(x=species,y=Structure,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Overall
p4 <-merge_gift_tm5 %>%
  ggplot(aes(x=species,y=Overall,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

Warning: Computation failed in `stat_compare_means()`.
Computation failed in `stat_compare_means()`.
Computation failed in `stat_compare_means()`.
Computation failed in `stat_compare_means()`.
Caused by error:
! argument "x" is missing, with no default

7.9.4 Comparison of samples from the 6th Time_point (6_Post-FMT2)

#Merge the functional domains with the metadata
merge_gift_TM6 <- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_TM6, by="Tube_code")

#Biosynthesis
p1 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Biosynthesis,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Type")

#Degradation
p2 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Degradation,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "type")

#Structure
p3 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Structure,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "type")

#Overall
p4 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Overall,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "type")

---

1. University of Copenhagen, garazi.bideguren@sund.ku.dk↩︎

2. University of Copenhagen, ostaizka.aizpurua@sund.ku.dk↩︎

3. Sociedad de Ciencias Aranzadi-Departamento de Herpetología, ccabido@aranzadi.eus↩︎

4. University of Copenhagen, antton.alberdi@sund.ku.dk↩︎